clc;
clear;

%% 迭代步骤 A 开始
% 读取网格数据
load msh

%% 迭代步骤 B 开始
% 设定物料黏度
niandu = 1000;

%% 迭代步骤 C 开始
% 设定边界条件
u1 = 0;
v1 = 0;
u3 = 0.01;
v3 = 0;
JBV1 = [BP1', u1*ones(size(BP1))', v1*ones(size(BP1))'];
JBV3 = [BP3', u3*ones(size(BP3))', v3*ones(size(BP3))'];
JBV = [JBV1; JBV3];
P2 = 0;
P4 = 1000;
JBP2 = [BE2, ones(size(BE2(:,1)))*P2, ones(size(BE2(:,1)))*P2];
JBP4 = [BE4, ones(size(BE4(:,1)))*P4, ones(size(BE4(:,1)))*P4];
JBP = [JBP2; JBP4];
clear JBV1 JBV3 BP1 BP3 BP4;
clear JBP2 JBP4 P2 P4;
clear BE1 BE2 BE3 BE4;
clear u1 v1 u3 v3;

%% 迭代步骤 D 开始
% 初始化总体方程各个数据子块
B1 = zeros(Nd, Nz);
B2 = zeros(Nd, Nz);
D11 = zeros(Nz, Nz);
D12 = zeros(Nz, Nz);
D21 = zeros(Nz, Nz);
D22 = zeros(Nz, Nz);
C1 = zeros(Nz, Nd);
C2 = zeros(Nz, Nd);
F1 = zeros(Nz, 1);
F2 = zeros(Nz, 1);

%% 迭代步骤 E 开始
% 计算单元方程系数矩阵各子块并组装
for i = 1:E
    for ie = 1:9
        JXYe(ie,:) = JXYV(JMV(i,ie),:);
    end
    [Be1,Be2] = function_of_Be(JXYe); % 调用 function_of_Be 函数
    [De11,De12,De21,De22] = function_of_De(JXYe, niandu); % 调用 function_of_De 函数
    [Ce1,Ce2] = function_of_Ce(JXYe); % 调用 function_of_Ce 函数
    for m = 1:4
        for n = 1:9
            B1(JMP(i,m),JMV(i,n)) = B1(JMP(i,m),JMV(i,n)) + Be1(m,n);
            B2(JMP(i,m),JMV(i,n)) = B2(JMP(i,m),JMV(i,n)) + Be2(m,n);
        end
    end
    for m = 1:9
        for n = 1:9
            D11(JMV(i,m),JMV(i,n)) = D11(JMV(i,m),JMV(i,n)) + De11(m,n);
            D12(JMV(i,m),JMV(i,n)) = D12(JMV(i,m),JMV(i,n)) + De12(m,n);
            D21(JMV(i,m),JMV(i,n)) = D21(JMV(i,m),JMV(i,n)) + De21(m,n);
            D22(JMV(i,m),JMV(i,n)) = D22(JMV(i,m),JMV(i,n)) + De22(m,n);
        end
    end
    for m = 1:9
        for n = 1:4
            C1(JMV(i,m),JMP(i,n)) = C1(JMV(i,m),JMP(i,n)) + Ce1(m,n);
            C2(JMV(i,m),JMP(i,n)) = C2(JMV(i,m),JMP(i,n)) + Ce2(m,n);
        end
    end
end

%% 迭代步骤 F 开始
% 计算单元方程右边向量子块并组装
for i = 1:length(JBP(:,1))
    for ie = 1:9
        JXYe(ie,:) = JXYV(JMV(JBP(i,1),ie),:);
    end
    [Fe1,Fe2] = function_of_Fe(JXYe, JBP(i,:)); % 调用 function_of_Fe 函数
    for m = 1:9
        F1(JMV(JBP(i,1),m),1) = F1(JMV(JBP(i,1),m),1) + Fe1(m,1);
        F2(JMV(JBP(i,1),m),1) = F2(JMV(JBP(i,1),m),1) + Fe2(m,1);
    end
end

%% 迭代步骤 G 开始
% 组合总体方程
K = [D11 D12 -C1;
    D21 D22 -C2;
    B1 B2 zeros(Nd,Nd)];
B = [-F1; -F2; zeros(Nd,1)];

%% 迭代步骤 H 开始
% 代入 JBV 数据
N_matrix = 2*Nz + Nd;
for i = 1:length(JBV(:,1))
    II = JBV(i,1);
    u = JBV(i,2);
    for J = 1:N_matrix
        B(J) = B(J) - K(J,II)*u;
    end
    K(II,:) = zeros(1,N_matrix);
    K(:,II) = zeros(N_matrix,1);
    K(II,II) = 1;
    B(II) = u;
end
for i = 1:length(JBV(:,1))
    II = Nz + JBV(i,1);
    v = JBV(i,3);
    for J = 1:N_matrix
        B(J) = B(J) - K(J,II)*v;
    end
    K(II,:) = zeros(1,N_matrix);
    K(:,II) = zeros(N_matrix,1);
    K(II,II) = 1;
    B(II) = v;
end

%% 迭代步骤 I 开始
% 清理内存, 求解方程
clear D11 D12 D21 D22 C1 C2 B1 B2
clear F1 F2 De11 De12 De21 De22 JXYe
clear Be1 Be2 Ce1 Ce2 JBP JMP JXYP
clear Fe1 Fe2 P_element P_size P_value
clear i i_JBP ie l_cos_theta_x m_cos_theta_y r s
x = K\B;
ux_k_1 = x(1:Nz);
vy_k_1 = x(1+Nz:2*Nz);
p4 = x(1+2*Nz:2*Nz+Nd);
p_k_1 = [Pding2Pzong(p4,JMV)]'; % 压力插值计算


%% 迭代步骤 J 开始
% 输出Tecplot后处理结果
E = E*4
Nz = Nz
data = [JXYV, ux_k_1, vy_k_1, sqrt(ux_k_1.^2 + vy_k_1.^2), p_k_1]
JMV4 = JMV_9to4(JMV)
% 清除多余变量
clear B E II J JBV JMV JMV4 JXYV K N_matrix Nd Nz
clear data niandu p4 u v x
% 存储结果
save result_of_n1

%% 出口速度提取
for i = 1:length(BP2)
    UB2(i,1) = ux_k_1(BP2(i),1);
end
UB2 = UB2

